#ifndef _CENTERLINE_THREAD_H_
#define _CENTERLINE_THREAD_H_

#include "JJSkeleton.h"

#include "Common.h"
#include "VolumeSource.h"
#include "FilteringThreadBase.h"
#include "itkConnectedThresholdImageFilter.h"
#include "itkFastMarchingImageFilter.h"
#include "itkThresholdSegmentationLevelSetImageFilter.h"
#include "itkBinaryThresholdImageFilter.h"
#include "itkCastImageFilter.h"
#include "itkImageFileWriter.h"
#include "itkCurvatureFlowImageFilter.h"
#include "itkShrinkImageFilter.h"
#include "itkExtractImageFilter.h"
#include "itkCurvatureAnisotropicDiffusionImageFilter.h"
#include "itkTimeProbe.h"
#include "itkResampleImageFilter.h"
#include "itkIdentityTransform.h"
#include "itkImageRegionIterator.h"
#include "itkNearestNeighborInterpolateImageFunction.h"
#include "itkConfidenceConnectedImageFilter.h"
#include "itkInvertIntensityImageFilter.h"
#include "itkBinaryContourImageFilter.h"
#include "itkBinaryPruningImageFilter.h"

#include <QList>
#include <QTest>



template<class CInputPixelType>
class CenterLineThread : public FilteringThreadBase
{	
	typedef itk::Image<CInputPixelType, 2> CInputImageType2D;
	typedef itk::Image<CInputPixelType, 3> CInputImageType;
	typedef itk::Image<float, 3> CImageTypeFloat3D;

	typedef itk::ImageDuplicator<CInputImageType> ImageDuplicator;
	typedef itk::ConnectedThresholdImageFilter<CInputImageType, CInputImageType> CTFilterType;
	typedef itk::ConfidenceConnectedImageFilter<CInputImageType, CInputImageType> CCFilterType;
	typedef itk::BinaryContourImageFilter<CInputImageType, CInputImageType> BinaryContourImageFilter;

public:
	CenterLineThread(Source *source):FilteringThreadBase(source)
	{
		m_replaceValue = 1;

		m_centerLineExtractor = new Skeleton();
		m_inputImage = ((VolumeSource<CInputImageType>*)m_source)->getVolume();

		m_biggestBranchOnly = false;
		m_pVesselEdge = NULL;
		m_pVesselVolume = NULL;
		m_pSkeleton = NULL;
		m_pControlPoints = NULL;
	}
	~CenterLineThread()
	{
	}

	void SetBiggestBranch(bool flag)
	{
		m_biggestBranchOnly = flag;
	}

	void Release()
	{
		if ( m_centerLineExtractor ) delete m_centerLineExtractor; m_centerLineExtractor = NULL;
		if ( m_pVesselEdge) delete m_pVesselEdge; m_pVesselEdge = NULL;
		if ( m_pVesselVolume) delete m_pVesselVolume; m_pVesselVolume = NULL;
		if ( m_pSkeleton) delete m_pSkeleton; m_pSkeleton = NULL;
		if ( m_pControlPoints) delete m_pControlPoints; m_pControlPoints = NULL;
	}

	virtual void addSeed(typename CInputImageType::IndexType seed)
	{
		m_seedList.push_back(seed);
	}

	virtual void setReplaceValue(int value = 1)
	{
		m_replaceValue = value;
	}

	int GetBranchAndExtremePoints(std::vector<Node> &vec)
	{
		return m_centerLineExtractor->GetBranchAndExtremePoints(vec);
	}

	typename CImageTypeFloat3D::Pointer GetDistanceMap()
	{
		return m_volume_DistanceMap;
	}

	void GetSpacing(double *spacing)
	{
		spacing[0] = this->spacing[0];
		spacing[1] = this->spacing[1];
		spacing[2] = this->spacing[2];
	}


protected:

	virtual void run()
	{
		typename CInputImageType::Pointer segmented;

		centerLineExtraction();
	}


	// combine with confidenceConnectedFiltering
	typename CInputImageType::Pointer centerLineExtraction()
	{
		ProgressCommandType::Pointer progressCmd = ProgressCommandType::New();
		progressCmd->SetCallbackFunction( this, &CenterLineThread::ProcessEvent );

		emit emitProgressStart(QString("Center Line Extraction..."));

		typename BinaryContourImageFilter::Pointer binContourFilter = BinaryContourImageFilter::New();
		binContourFilter->SetInput(m_inputImage);
		binContourFilter->SetForegroundValue(1);
		binContourFilter->SetBackgroundValue(0);
		binContourFilter->Update();

		if ( m_pVesselEdge) delete m_pVesselEdge;
		m_pVesselEdge = new RxMaskVolume();
		if ( m_pVesselVolume) delete m_pVesselVolume;
		m_pVesselVolume = new RxMaskVolume();
		if ( m_pSkeleton) delete m_pSkeleton;
		m_pSkeleton = new RxMaskVolume();
		if ( m_pControlPoints) delete m_pControlPoints;
		m_pControlPoints = new RxMaskVolume();

		typename CInputImageType::SizeType size = m_inputImage->GetBufferedRegion().GetSize();
		typename CInputImageType::SpacingType spac = m_inputImage->GetSpacing();

		m_pVesselEdge->CreateMaskVolume(size[0]+4, size[1]+4, size[2]+4);
		m_pVesselVolume->CreateMaskVolume(size[0]+4, size[1]+4, size[2]+4);
		m_pSkeleton->CreateMaskVolume(size[0]+4, size[1]+4, size[2]+4);
		m_pControlPoints->CreateMaskVolume(size[0]+4, size[1]+4, size[2]+4);

		spacing[0] = spac[0];
		spacing[1] = spac[1];
		spacing[2] = spac[2];

		typename CInputImageType::IndexType pixelIndex;
		for (unsigned int i = 0; i < size[2]; i++)
		{
			for (unsigned int j = 0; j< size[1]; j++ )
			{
				for ( unsigned int k = 0; k < size[0]; k++)
				{
					pixelIndex[0] = k;
					pixelIndex[1] = j;
					pixelIndex[2] = i;
					if ( m_inputImage->GetPixel(pixelIndex) > 0)
						m_pVesselVolume->SetMask(k+2, j+2, i+2);

					if ( binContourFilter->GetOutput()->GetPixel(pixelIndex) == 1)
					{
						m_pVesselEdge->SetMask(k+2, j+2, i+2);
					}
				}
			}
		}


 		m_centerLineExtractor->SetParameter(m_pVesselEdge, m_pVesselVolume, m_pSkeleton, m_pControlPoints, 
			0, size[0]-1, 0, size[1]-1, 0, size[2]-1, false);
		m_centerLineExtractor->SetBiggestBranch(m_biggestBranchOnly);
 		m_centerLineExtractor->Skeletonization();
		//float *pDistanceMapFromPS = m_centerLineExtractor->GetDistanceMapFromPS();

		typename ImageDuplicator::Pointer dup_VesselEdge = ImageDuplicator::New();
		typename ImageDuplicator::Pointer dup_VesselVolume = ImageDuplicator::New();
		typename ImageDuplicator::Pointer dup_Skeleton = ImageDuplicator::New();
		typename ImageDuplicator::Pointer dup_ControlPoints = ImageDuplicator::New();
		//ImageDuplicator::Pointer dup_DistanceMap = ImageDuplicator::New();

		dup_VesselEdge->SetInputImage(m_inputImage);dup_VesselEdge->Update();
		dup_VesselVolume->SetInputImage(m_inputImage);dup_VesselVolume->Update();
		dup_Skeleton->SetInputImage(m_inputImage);dup_Skeleton->Update();
		dup_ControlPoints->SetInputImage(m_inputImage);	dup_ControlPoints->Update();
		//dup_DistanceMap->SetInputImage(m_inputImage);	dup_DistanceMap->Update();

		typename CInputImageType::Pointer volume_VesselEdge = dup_VesselEdge->GetOutput(); volume_VesselEdge->FillBuffer(0);
		typename CInputImageType::Pointer volume_VesselVolume = dup_VesselVolume->GetOutput(); volume_VesselVolume->FillBuffer(0);
		typename CInputImageType::Pointer volume_Skeleton = dup_Skeleton->GetOutput(); volume_Skeleton->FillBuffer(0);
		typename CInputImageType::Pointer volume_ControlPoints = dup_ControlPoints->GetOutput(); volume_ControlPoints->FillBuffer(0);
		//m_volume_DistanceMap = dup_DistanceMap->GetOutput(); volume_DistanceMap->FillBuffer(0);

		m_volume_DistanceMap = CImageTypeFloat3D::New();
		m_volume_DistanceMap->SetRegions(m_inputImage->GetLargestPossibleRegion());
		m_volume_DistanceMap->Allocate();
		m_volume_DistanceMap->FillBuffer(0);

		for (unsigned int i = 0; i < size[2]; i++)
		{
			for (unsigned int j = 0; j< size[1]; j++ )
				for ( unsigned int k = 0; k < size[0]; k++)
				{
					pixelIndex[0] = k;
					pixelIndex[1] = j;
					pixelIndex[2] = i;

					if (m_pVesselEdge->GetMask(k+2,j+2,i+2) != 0)
						volume_VesselEdge->SetPixel(pixelIndex, 1);

					if (m_pSkeleton->GetMask(k+2,j+2,i+2) != 0)
						volume_Skeleton->SetPixel(pixelIndex, 1);

					if (m_pControlPoints->GetMask(k+2,j+2,i+2) != 0)
						volume_ControlPoints->SetPixel(pixelIndex, 1);

					m_volume_DistanceMap->SetPixel(pixelIndex, m_centerLineExtractor->GetDistanceMapValue(k+2, j+2, i+2));
				}
		}

// 		VolumeSource<CInputImageType> *subsource2 = 
// 			new VolumeSource<CInputImageType>(QString("ExtremePoints#%1").arg(m_source->rowCount()+1));
// 		subsource2->setVolume(volume_VesselEdge);
// 		m_source->addChild(subsource2);


		VolumeSource<CInputImageType> *subsource4 = 
			new VolumeSource<CInputImageType>(QString("Skeleton#%1").arg(m_source->rowCount()+1));
		subsource4->setVolume(volume_Skeleton);
		m_source->addChild(subsource4);


// 		VolumeSource<CInputImageType> *subsource5 = 
// 			new VolumeSource<CInputImageType>(QString("ControlPoints#%1").arg(m_source->rowCount()+1));
// 		subsource5->setVolume(volume_ControlPoints);
// 		m_source->addChild(subsource5);
// 
// 		VolumeSource<CImageTypeFloat3D> *m_distanceMapFromPS = new VolumeSource<CImageTypeFloat3D>(QString("DistanceMap"));
// 		m_distanceMapFromPS->setVolume(m_volume_DistanceMap);
// 		m_source->addChild(m_distanceMapFromPS);
// 		filter->AddObserver( itk::ProgressEvent(), progressCmd.GetPointer());
// 		filter->Update();

		emit emitProgressDone();

		return m_inputImage;
	}

protected:

	int m_replaceValue;
	double spacing[3];
	bool m_biggestBranchOnly;
	
	Skeleton *m_centerLineExtractor;

	QList<typename CInputImageType::IndexType> m_seedList;
	typename CInputImageType::Pointer m_inputImage;
	typename CImageTypeFloat3D::Pointer m_volume_DistanceMap;

	RxMaskVolume *m_pVesselEdge;
	RxMaskVolume *m_pVesselVolume;
	RxMaskVolume *m_pSkeleton;
	RxMaskVolume *m_pControlPoints;
};






#endif
